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EXECUTIVE  SUMMARY 


The  time  series  representing  the  acoustic  energy  received  at  two  sensors  due  to  a  source  moving 
through  an  oceanic  waveguide  can  exhibit  strong  spectral  correlation  between  different  frequencies. 
This  correlation  can  be  exploited  to  enhance  passive  detection  of  moving  targets — particularly  in 
shallow  water  and  near-field  environments. 

The  structure  of  this  spectral  correlation  is  derived  in  this  document.  It  is  also  shown  that  estab¬ 
lished  techniques  can  estimate  the  spectral  correlation  from  time  series.  This  algorithm  is  assessed  on 
real  underwater  acoustic  data.  Excellent  agreement  is  observed  with  theory. 
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1.  INTRODUCTION 


Many  algorithms  used  in  SONAR  signal  processing  require  the  estimation  of  auto-  and  cross¬ 
spectra  of  time  series  received  at  an  array  of  sensors.  A  fundamental  assumption  is  that  the  time 
series  is  stationary.  However,  a  moving  source  destroys  stationarity  between  sensors.  This  is  quite 
problematic  for  the  estimation  of  second-order  moments.  In  practice,  one  assumes  that  the  data  are 
piecewise  stationary,  i.e.,  that  the  moments  can  be  estimated  from  sufficiently  long  segments  of  data. 
This  assumption  may  not  be  valid  for  sources  in  a  near-field  or  shallow-water  environment. 

In  this  document,  we  define  and  derive  the  2-D  spectrum  of  a  broadband  source  moving  through 
an  oceanic  waveguide.  We  show  how  an  established  technique  exploits  the  structure  of  the  2-D  spec¬ 
trum  to  estimate  a  related  quantity,  the  2-D  spectral  coherence.  This  technique  accounts  for  source 
motion  in  the  spectral  analysis  algorithm  by  using  the  spectral  correlation  properties  of  non-station- 
ary  processes.  By  incorporating  this  spectral  analysis  approach,  the  detection  performance  of  signal 
and  array  processing  algorithms  is  enhanced  since  the  2-D  spectral  coherence  can  be  used  as  a  dis¬ 
criminant  for  moving  targets. 

In  section  2,  we  derive  the  2-D  spectrum  for  an  acoustic  source  moving  through  an  oceanic  wave¬ 
guide  and  describe  the  algorithm  that  estimates  the  2-D  spectral  coherence  from  time  series.  In  sec¬ 
tion  3,  we  move  from  this  theory  to  the  real  world.  We  analyze  experimental  data  and  compare  these 
results  to  the  theory.  Section  4  contains  a  summary  and  in  section  5  we  recommend  several  ways  to 
exploit  these  results. 


2.  2-D  SPECTRUM  OF  A  MOVING  SOURCE 

Hawker  (1979)  originally  derived  the  pressure  field  generated  by  an  acoustic  source  moving 
through  an  oceanic  waveguide.  We  related  Hawker’s  results  to  time-varying  linear  system  theory  and 
used  them  as  a  basis  for  an  acoustic  time  series  simulation  (Reuter,  1994).  In  this  section,  we  briefly 
review  these  results.  We  then  derive  the  2-D  spectrum  (Hurd,  1989)  of  an  acoustic  source  moving 
through  a  waveguide  and  describe  an  algorithm  that  generates  a  consistent  estimate  of  the  2-D  spec¬ 
tral  coherence  (Hurd  &  Gerr,  1991). 

2.1  2-D  SPECTRUM  IN  TIME-VARYING  LINEAR  SYSTEMS 

The  impulse  response  and  the  bifrequency  system  function  r{o),  u)  of  a  time-varying  linear  system 
are  related  by  the  Fourier  transform  pair  (Reuter,  1994;  Zadeh,  1950) 

h{t,  —  t)  r{(jD,  u).  (1) 

If  {yi(0}  and  {y2(0}  are  the  outputs  of  two  real,  stable  linear  systems  driven  by  the  same  real,  wide 
sense  stationary  process  {x(0) ,  then  the  2-D  spectrum  between  {yi(0}  and  {y2(0}  is  given  by  Reuter 
(1994) 


^  j"  SJfi))r^{a),a)r*2{a},p)d(o, 


(2) 


1 


where  Sx(o))  is  the  power  spectral  density  of  the  process  {x(0)-  5'j2(ct,y3)  describes  the  joint  spectral 
correlation  property  of  processes  {yi(0}  and  {y2(0}  • 

If  the  linear  systems  were  also  time  invariant,  the  2-D  spectrum  can  be  shown  to  be  (Reuter,  1994) 

5i2(a,y8)  =  2jtSx0)H^0)H*2(fi)d0  -  a),  (3) 

where 

H{a)\t)  =  j  r(a),u)e““du 

and  d{a)  is  the  Dirac  delta  function.  Since  the  processes  {yi(0}  and  {y2(0}  are  jointly  wide  sense  sta¬ 
tionary,  there  is  no  correlation  between  different  frequencies.  This  is  also  demonstrated  in  equa¬ 
tion  (3)  where  the  support  of  the  2-D  spectrum  5'i2(«,^)  is  only  on  the  a  =  fi  portion  of  the  spectral 
coherence  plane. 


2.2  2-D  SPECTRUM  OF  A  SOURCE  MOVING  THROUGH  AN  OCEANIC  WAVEGUIDE 

In  order  to  derive  a  useful  bifrequency  system  function  of  an  acoustic  source  moving  through  an 
oceanic  waveguide  with  respect  to  a  sensor.  Hawker  makes  the  linear  motion  model  approximation. 
The  schematic  of  this  system  is  shown  in  figure  1 ,  where  R{t)  and  6{t)  are  time-dependent  range  and 
angles.  In  this  approximation,  the  source  is  moving  sufficiently  slowly  so  that  R{t)  and  sin  {6{t))  can 
be  adequately  represented  by  the  linear  terms  of  the  expansions  about  some  region  of  time,  A  t,  cen¬ 
tered  at  t  —  to  as 


R(t)  ==  R(to)  +  (t-to)R'(to) 

and  (5) 

sin  ~  sin  (6{to))  +  {t-to)d'{to)  cos  {0{to). 

Using  the  assumption  that  the  acoustic  propagation  environment  is  cylindrically  symmetric. 
Hawker  shows  that  the  bifrequency  system  function  is  given  by 

r((o,u)  = 


ik„RA  1  -  ;j^jSin  (d„) 


(v„Ay2)e 


i{Aii  —  (D)M 


D 


(6) 


where 

Zn  =  n*^  modal  depth  eigenfunction 
kn  =  n^^  modal  eigenvalue 
=  fith  modal  group  velocity 
an  =  modal  attenuation  term 
V  =  source  speed 
Zs  =  source  depth 

Mj)  =  delay  associated  with  first  modal  arrival  time. 
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Figure  1 .  Source/receiver  track  scenario. 


with  Ro  —  R{to) ,  9o  —  0{to) ,  An  =  knV  sin  {6o)  —  and  =  m  —  a>  +  A„.  For  notational 

simplicity,  the  explicit  dependence  of  these  variables  on  o)  is  suppressed. 

First,  without  loss  of  generality,  we  set  to  —  0.  For  analysis  purposes,  we  let  Zl  ^  ^  oo  in  equa¬ 
tion  (6).  Then,  using  the  identity 

^  00  sinc(v„  Af/2)  =  27td{Vn) ,  * 


we  get 


r{a),u)=  exp 

V  2  y 


-iknRo  I  1  -  sin  {Oo) 


CT(V  le 


Now,  using  equation  (8)  in  equation  (2)  and  grouping  components  dependent  on  cd  into  a  modal 
dependent  term  En^^iiwy  we  get 


512(0:, /?) 


n  m 


En,m{0))d{a  -  0)  +  A^n)^{fi  -  (o  +  A2jdo), 


where  the  modal  sums  in  equation  (9)  are  assumed  to  be  finite.  Next,  let 


(9) 


fn,m(0:)/5) 


j  En,„r(0))d(a  -  0)  +  A^n)^(^  -  0)  +  A2jdQ). 


(10) 


Proceeding  as  in  Papoulis  (1962),  we  interpret  In,m{o.,f^)  as  a  distribution.  Then  we  evaluate  the  inte¬ 
gral 


I  In,m{a,p)(t>i^)d^ 


En^J((i))d{a  —  (o  +  A^n)^0  -  0)  +  A2„)(t)(l3)dcDd(i, 


(11) 
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where  (j){fi)  is  a  continuous  test  function.  After  exchanging  the  order  of  integration,  we  get 


=  E„^^(co)d{a  -  (0  +  Ai„)0((u  -  A2jd(o. 


(12) 


Next,  we  approximate  and  A  2^  with  the  linear  terms  of  the  Taylor  series  expansion  around  a 
frequency  (Oc  within  the  band  of  interest  as 


^  +  A\„((Oc)(a)  -  COc). 

Then,  using  the  fact  that  A^^((Oc)  —  A'i„(cUc)a>c  =  0  and  after  a  little  manipulation  we  get 


1 


a 


1  -  a; 


2m 


1-A)J|  "ni-Ai 


[n  > 


1  -  a; 


a  \. 


In 


So,  from  equation  (14)  we  conclude  that 


fn,m(tt,/5)  ^  Efiifil 

I  1  “  ^In  II  \  1  “  ^In 


«  _U|  ^  ^2m 


1  -  A'l, 


\a  , 


where 


/  \  ^ 

A'in  =  ^gSiniOJ  +  I  (knvj«-  1). 

^  n  n  / 


And  so  the  2-D  spectrum  is  given  by 


(14) 


(15) 


(16) 


(17) 


We  see  from  equations  (15)  and  (17)  that  Si2(c£,yS)  only  has  support  on  regions  of  the  (a,^)- 

spectral  plane  where  =  (1  —  A2m)/(1  ~  A'i„)a.  From  equation  (16)  we  see  that  the  slopes  of  these 
spectral  lines  are  dependent  on  the  source  velocity  with  respect  to  each  sensor  as  well  as  the  acoustic 
propagation  environment.  The  cross-modal  spectral  correlation  is  also  evident. 


2.3  2-D  SPECTRAL  COHERENCE 

The  structure  of  the  2-D  spectrum  in  equations  (15)  and  (17)  lends  itself  to  an  estimation  proce¬ 
dure  first  introduced  by  Hurd  &  Gerr  (1991)  for  the  spectral  analysis  of  cyclostationary  processes 
and  adapted  by  Allen  to  include  the  Doppler  processing  of  moving  targets  in  a  plane-wave  propaga¬ 
tion  environment  (Allen  &  Hobbs,  1992,  1993). 

The  procedure  involves  forming  the  outer  product  of  the  short-time  Fourier  transform  of  the 
sampled  data  from  each  sensor  and  smoothing  along  lines  corresponding  to  potential  source-sensor 
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trajectories.  This  method  forms  a  consistent  estimate  5j2(a,y3)  of  the  sum  of  the  coefficients  of  the 
delta  functions  in  equation  (15).  The  2-D  spectral  coherence  then  is  defined  as  (Gerr  &  Allen,  1994). 


Cj2(ct,/1) 


a)5'22(/?,^) 


(18) 


3.  EXPERIMENTAL  RESULTS 

In  this  section  we  demonstrate  the  spectral  analysis  technique  from  section  2.3  on  a  real  data  set. 
We  analyzed  a  portion  of  data  collected  during  the  SWellEX-3  experiment  that  was  held  during  July 
and  August  of  1994  in  the  Southern  California  Shallow-Water  Testbed  (SWTB).  This  experiment 
was  meant  to  test  the  passive  detection  of  targets  in  a  shallow-water  environment  (50  to  200  m). 

3.1  DATA  SET  DESCRIPTION 

The  SWellEX-3  data  that  were  analyzed  were  part  of  the  Short-Range  CPA  Event.  We  used  30 
minutes  of  data  collected  from  hydrophone  2  of  node  2  (sensor  1)  and  hydrophone  131  of  node  5 
(sensor  2)  of  the  Shallow-Water  Sensor  String  (SWSS)  horizontal  line  array.  Both  hydrophones  are 
of  type  HE(5-1040  Hz}).  The  bottom-mounted  array  was  deployed  at  a  depth  of  approximately 
170  to  180  m.  The  environmental  information  is  displayed  in  appendix  A  (Baxley  &  Booth,  1995)  in 
the  form  of  a  KRAKEN  (Porter,  1991)  environmental  input  file  for  an  adjacent  region.  Eigure  2  is  a 
schematic  of  the  data  set,  which  includes  the  relative  sensor  locations  and  the  source  track  deter¬ 
mined  from  Global  Positioning  System  (GPS)  data.  The  broadband  source  was  produced  by  a  pseu¬ 
do-random  noise  generator  and  was  towed  at  a  speed  of  approximately  1.54  m/s  at  a  nominal  depth 
of  60  m.  The  spectrograms  of  the  two  hydrophones  are  given  in  figures  3  and  4  over  the  400-  to 
800-Hz  frequency  range. 

3.2  SPECTRAL  ANALYSIS  RESULTS 

We  present  results  from  three  principle  subsections  of  this  data  set  labeled  A,  B,  and  C  in  figure  2. 
They  represent  varying  degrees  of  relative  range-rate  with  respect  to  the  sensors.  Section  A  is  at  the 
beginning  of  the  data  set  and  corresponds  to  the  least  amount  of  relative  Doppler  offset  between  the 
two  sensors.  Section  C  represents  the  most  severe  case.  Appendix  B  contains  the  MATEAB  code 
used  to  calculate  the  2-D  spectral  coherence. 

In  the  spectral  analysis  results  presented  in  this  section,  it  was  necessary  to  eliminate  an  estimate 
of  the  “coarse  time  delay”  associated  with  the  difference  in  propagation  time  between  the  two  sen¬ 
sors.  This  is  because  this  estimation  procedure  is  based  upon  a  spectral  smoothing  approach  and  so  is 
very  sensitive  to  the  phase  introduced  by  the  propagation  delay  (Gerr  &  Allen,  1994). 

Eigures  5,7,  and  9  are  the  2-D  spectral  coherence  for  the  three  data  sections  over  the  660-  to  670- 
Hz  band  generated  from  15.7  seconds  of  time  series.  The  a-coordinate  is  on  the  vertical  axis.  To 
more  clearly  display  the  frequency  offset  of  the  spectral  correlation  structure,  figures  11  and  13  dis¬ 
play  the  2-D  spectral  coherence  of  data  sections  B  and  C  over  the  660-  to  665-Hz  frequency  band. 
Eigures  6,  8,  10,  12,  and  14  are  the  corresponding  analytical  regions  of  support  generated  with  the 
results  of  section  2  using  22  modes  and  using  KRAKEN  with  the  environmental  information  given  in 
appendix  A.  Only  the  potential  regions  of  support  are  plotted  using  equation  (15)  and  not  the  actual 
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Figure  2.  Source  track  and  sensor  geometry. 
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Figure  9.  2-D  spectral  coherence  for  section  C  over 
660-  to  670-Hz  band. 


Figure  10.  Analytical  region  of  support  of  2-D  spectral 
coherence  for  section  C  over  660-  to  670-Hz  band. 
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Figure  11.  2-D  spectral  coherence  for  section  B 
660-  to  665-Hz  band. 


Figure  12.  Analytical  region  of  support  of  2-D  s 
coherence  for  section  B  over  660-  to  665-Hz  ba 
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Figure  13.  2-D  spectral  coherence  for  section  C  over 
660-  to  665-Hz  band. 
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Figure  14.  Analytical  region  of  support  of  2-D  spectral 
coherence  for  section  C  over  660-  to  665-Hz  band. 
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2-D  spectral  coherence.  Since  no  summing  is  being  done,  many  more  potential  cross-modal  lines  are 
plotted  in  these  figures.  However,  they  are  in  approximate  agreement  with  the  corresponding  esti¬ 
mated  plots.  Discrepancies  are  probably  due  to  the  inaccuracy  of  the  environmental  information  and 
the  fact  that  the  actual  propagation  environment  for  this  site  is  not  range-independent  as  assumed  in 
the  analysis  of  section  2.  The  spectral  correlation  due  to  source  motion  is  clearly  evident  in  the  plots. 

One  of  the  primary  attributes  of  spectral  coherence  processing  is  that  it  allows  for  very  long  Four¬ 
ier  transforms  to  be  taken.  This  is  in  contrast  to  conventional  spectral  estimation  that  requires  the 
source  to  be  approximately  stationary  over  the  data  length  of  the  transform.  Since  source  motion  is 
inherently  accounted  for  in  the  spectral  coherence  estimation  procedure,  this  type  of  processing 
accommodates  much  longer  transform  lengths.  Figure  15  is  the  estimated  2-D  spectral  coherence 
over  the  660-  to  660.625-Hz  frequency  band  for  data  section  A  generated  from  251.7  seconds  of  time 
series.  Figure  16  is  the  corresponding  analytical  region  of  support.  Much  finer  spectral  coherence 
detail  is  evident  in  these  plots  as  opposed  to  the  detail  afforded  by  the  15.7-second  transform  repre¬ 
sented  in  figures  5  and  6  over  the  wider  frequency  range.  We  see  interesting  auto-  and  cross-modal 
spectral  correlation  from  these  plots.  However,  this  length  of  Fourier  transform  would  be  inappropri¬ 
ate  for  the  shorter-range  sections,  B  and  C,  since  the  linear  motion  model  would  be  invalid  in  these 
sections  over  this  length  of  time. 
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4.  CONCLUSION 


In  this  document  we  derived  the  2-D  spectrum  of  a  broadband  source  moving  through  an  oceanic 
waveguide.  The  2-D  spectrum  describes  the  spectral  correlation  structure  of  the  time  series  received 
at  two  sensors  due  to  a  moving  source.  We  demonstrated  how  its  structure  is  dependent  on  the  source 
velocity  with  respect  to  the  sensors  as  well  as  the  acoustic  propagation  environment. 

We  also  showed  that  this  structure  lends  itself  to  an  established  spectral  smoothing  procedure  that 
was  shown  to  be  a  consistent  estimator  of  the  2-D  spectral  coherence.  We  then  applied  this  technique 
to  a  real  data  set  collected  during  the  SWellEX-3  experiment  and  demonstrated  how  the  theoretical 
spectral  coherence  approximates  the  experimental. 


5.  RECOMMENDATIONS 

The  most  critical  issue  in  using  the  results  presented  in  this  document  is  the  estimation  procedure. 
Further  work  should  be  done  to  test  the  robustness  of  the  spectral  smoothing  technique.  For  instance, 
it  should  be  tested  in  an  environment  in  which  there  is  a  strong  interference  overlapping  the  band  of 
the  desired  source.  The  effect  of  near-  and  far-field  interferers  on  the  estimation  technique  should  be 
determined. 

Several  improvements  can  be  made  to  the  procedure  applied  to  the  underwater  acoustic  problem. 
An  automated  procedure  should  be  developed  that  would  accommodate  many  source  ranges  and 
velocities  and  simultaneously  account  for  coarse  time  delay  correction.  Also,  some  effort  should  be 
expended  to  develop  a  computationally  efficient  implementation. 

Finally,  a  method  can  be  developed  to  incorporate  2-D  spectral  coherence  into  an  array-processing 
algorithm  to  enhance  the  detection  and  localization  of  moving  targets.  This  should  be  completely 
analogous  to  power  spectral  density  use  in  frequency-domain,  array-processing  algorithms.  The 
theoretical  2-D  spectrum  should  also  be  derived  for  a  range-dependent  propagation  environment. 
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APPENDIX  A 

EXPERIMENT  SITE  ENVIRONMENTAL  INFORMATION 


This  appendix  contains  the  nominal  KRAKEN  environmental  file  for  the  SWellEX-3  experiment. 


"SHELLEX-3" 

10 

3 

'NVF.  ' 

400  0.0  198.0 


0.0 

1522.610 

0.0 

1.0 

4.55578E-06 

0.0 

47.0 

1497.400 

0.0 

1.0 

4.55578E-06 

0.0 

1.0 

1523.110 

0.0 

1.0 

4.55578E-06 

0.0 

48.0 

1497.150 

0.0 

1.0 

4.55578E-06 

0.0 

2.0 

1523.180 

0.0 

1.0 

4.55578E-06 

0.0 

49.0 

1497.050 

0.0 

1.0 

4.55578E-06 

0.0 

3.0 

1523.230 

0.0 

1.0 

4.55578E-06 

0.0 

50.0 

1496.860 

0.0 

1.0 

4.55578E-06 

0.0 

4.0 

1523.020 

0.0 

1.0 

4.55578E-06 

0.0 

51.0 

1496.360 

0.0 

1.0 

4.55578E-06 

0.0 

5.0 

1522.770 

0.0 

1.0 

4.55578E-06 

0.0 

52.0 

1496.280 

0.0 

1.0 

4.55578E-06 

0.0 

6.0 

1522.230 

0.0 

1.0 

4.55578E-06 

0.0 

53.0 

1496.170 

0.0 

1.0 

4.55578E-06 

0.0 

7.0 

1521.110 

0.0 

1.0 

4.55578E-06 

0.0 

54.0 

1495.960 

0.0 

1.0 

4.55578E-06 

0.0 

8.0 

1520.380 

0.0 

1.0 

4.55578E-06 

0.0 

55.0 

1495.600 

0.0 

1.0 

4.55578E-06 

0.0 

9.0 

1519.650 

0.0 

1.0 

4.55578E-06 

0.0 

56.0 

1495.610 

0.0 

1.0 

4.55578E-06 

0.0 

10.0 

1518.720 

0.0 

1.0 

4.55578E-06 

0.0 

57.0 

1495.550 

0.0 

1.0 

4.55578E-06 

0.0 

11.0 

1518.190 

0.0 

1.0 

4.55578E-06 

0.0 

58.0 

1495.330 

0.0 

1.0 

4.55578E-06 

0.0 

12.0 

1518.170 

0.0 

1.0 

4.55578E-06 

0.0 

59.0 

1494.940 

0.0 

1.0 

4.55578E-06 

0.0 

13.0 

1515.900 

0.0 

1.0 

4.55578E-06 

0.0 

60.0 

1494.690 

0.0 

1.0 

4.55578E-06 

0.0 

14.0 

1515.460 

0.0 

1.0 

4.55578E-06 

0.0 

61.0 

1494.550 

0.0 

1.0 

4.55578E-06 

0.0 

15.0 

1514. 190 

0.0 

1.0 

4.55578E-06 

0.0 

62.0 

1494.330 

0.0 

1.0 

4.55578E-06 

0.0 

16.0 

1514.040 

0.0 

1.0 

4.55578E-06 

0.0 

63.0 

1494.080 

0.0 

1.0 

4.55578E-06 

0.0 

17.0 

1513.300 

0.0 

1.0 

4.55578E-06 

0.0 

64.0 

1493.950 

0.0 

1.0 

4.55578E-06 

0.0 

18.0 

1511.930 

0.0 

1.0 

4.55S78E-06 

0.0 

65.0 

1494.160 

0.0 

1.0 

4.55578E-06 

0.0 

19.0 

1510.970 

0.0 

1.0 

4.55578E-06 

0.0 

66.0 

1493.880 

0.0 

1.0 

4.55578E-06 

0.0 

20.0 

1509.710 

0.0 

1.0 

4.55578E-06 

0.0 

67.0 

1494.100 

0.0 

1.0 

4.55578E-06 

0.0 

21.0 

1509.290 

0.0 

1.0 

4.55578E-06 

0.0 

68.0 

1493.970 

0.0 

1.0 

4.55578E-06 

0.0 

22.0 

1508.800 

0.0 

1.0 

4.55578E-06 

0.0 

69.0 

1493.880 

0.0 

1.0 

4.55578E-06 

0.0 

23.0 

1508.110 

0.0 

1.0 

4.55578E-06 

0.0 

70.0 

1493.940 

0.0 

1.0 

4.55578E-06 

0.0 

24.0 

1504.960 

0.0 

1.0 

4.55578E-06 

0.0 

71.0 

1493.580 

0.0 

1.0 

4.55578E-06 

0.0 

25.0 

1505.510 

0.0 

1.0 

4.55578E-06 

0.0 

72.0 

1493.420 

0.0 

1.0 

4.55578E-06 

0.0 

26.0 

1503.730 

0.0 

1.0 

4.55578E-06 

0.0 

73.0 

1494.060 

0.0 

1.0 

4.55578E-06 

0.0 

27.0 

1502.710 

0.0 

1.0 

4.55578E-06 

0.0 

74.0 

1493.600 

0.0 

1.0 

4.55578E-06 

0.0 

28.0 

1502.310 

0.0 

1.0 

4.55578E-06 

0.0 

75.0 

1493.550 

0.0 

1.0 

4.55578E-06 

0.0 

29.0 

1501.100 

0.0 

1.0 

4.55578E-06 

0.0 

76.0 

1493.050 

0.0 

1.0 

4.55578E-06 

0.0 

30.0 

1500.210 

0.0 

1.0 

4.55578E-06 

0.0 

77.0 

1493.010 

0.0 

1.0 

4.55578E-06 

0.0 

31.0 

1500.120 

0.0 

1.0 

4.55578E-06 

0.0 

78.0 

1492.780 

0.0 

1.0 

4.55578E-06 

0.0 

32.0 

1500.200 

0.0 

1.0 

4.55578E-06 

0.0 

79.0 

1492.560 

0.0 

1.0 

4.55578E-06 

0.0 

33.0 

1500.000 

0.0 

1.0 

4.55578E-06 

0.0 

80.0 

1492.280 

0.0 

1.0 

4.55578E-06 

0.0 

34.0 

1499.580 

0.0 

1.0 

4.55578E-06 

0.0 

81.0 

1492.210 

0.0 

1.0 

4.55578E-06 

0.0 
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94.0 

1491.470 

0.0 

1.0 

4.55578E-06 

0.0 

149.0 

1490.990 

0.0 

1.0 

4.55578E-06 

0.0 

95.0 

1491.280 

0.0 

1.0 

4.55578E-06 

0.0 

150.0 

1490.970 

0.0 

1.0 

4.55578E-06 

0.0 

96.0 

1491.000 

0.0 

1.0 

4.55578E-06 

0.0 

151.0 

1490.960 

0.0 

1.0 

4.55578E-06 

0.0 

97.0 

1490.870 

0.0 

1.0 

4.55578E-06 

0.0 

152.0 

1490.960 

0.0 

1.0 

4.55578E-06 

0.0 

98.0 

1490.930 

0.0 

1.0 

4.55578E-06 

0.0 

153.0 

1490.970 

0.0 

1.0 

4.55578E-06 

0.0 

99.0 

1491.060 

0.0 

1.0 

4.55578E-06 

0.0 

154.0 

1490.980 

0.0 

1.0 

4.55578E-06 

0.0 

100.0 

1491.260 

0.0 

1.0 

4.55578E-06 

0.0 

155.0 

1490.940 

0.0 

1.0 

4.55578E-06 

0.0 

101.0 

1491.180 

0.0 

1.0 

4.55578E-06 

0.0 

156.0 

1490.970 

0.0 

1.0 

4.55578E-06: 

0.0 

102.0 

1491.230 

0.0 

1.0 

4.55578E-06 

0.0 

157.0 

1490.980 

0.0 

1.0 

4.55578E-06 

0.0 

103.0 

1491.410 

0.0 

1.0 

4.55578E-06 

0.0 

158.0 

1490.960 

0.0 

1.0 

4.55578E-06 

0.0 

104.0 

1491.340 

0.0 

1.0 

4.55578E-06 

0.0 

159.0 

1490.940 

0.0 

1.0 

4.55578E-06 

0.0 

105.0 

1491.450 

0.0 

1.0 

4.55578E-06 

0.0 

160.0 

1490.980 

0.0 

1.0 

4.55578E-06 

0.0 

106.0 

1491.430 

0.0 

1.0 

4.55578E-06 

0.0 

161.0 

1490.980 

0.0 

1.0 

4.55578E-06 

0.0 

107.0 

1491.380 

0.0 

1.0 

4.55578E-06 

0.0 

162.0 

1490.950 

0.0 

1.0 

4.55578E-06 

0.0 

108.0 

1491.350 

0.0 

1.0 

4.55578E-06 

0.0 

163.0 

1490.960 

0.0 

1.0 

4.55578E-06 

0.0 

109.0 

1491.340 

0.0 

1.0 

4.55578E-06 

0.0 

164.0 

1490.940 

0.0 

1.0 

4.55578E-06 

0.0 

110.0 

1491.450 

0.0 

1.0 

4.55578E-06 

0.0 

165.0 

1490.940 

0.0 

1.0 

4.55578E-06 

0.0 

111.0 

1491.520 

0.0 

1.0 

4.55578E-06 

0.0 

166.0 

1490.950 

0.0 

1.0 

4.55578E-06 

0.0 

112.0 

1491.510 

0.0 

1.0 

4.55578E-06 

0.0 

167.0 

1490.950 

0.0 

1.0 

4.55578E-06 

0.0 

113.0 

1491.500 

0.0 

1.0 

4.55578E-06 

0.0 

168.0 

1490.870 

0.0 

1.0 

4.55578E-06 

0.0 

114.0 

1491.370 

0.0 

1.0 

4.55578E-06 

0.0 

169.0 

1490.860 

0.0 

1.0 

4.55578E-06 

0.0 

115.0 

1491.290 

0.0 

1.0 

4.55578E-06 

0.0 

170.0 

1490.860 

0.0 

1.0 

4.55578E-06 

0.0 

116.0 

1491.400 

0.0 

1.0 

4.55578E-06 

0.0 

171.0 

1490.840 

0.0 

1.0 

4.55578E-06 

0.0 

117.0 

1491.290 

0.0 

1.0 

4.55578E-06 

0.0 

172.0 

1490.810 

0.0 

1.0 

4.55578E-06 

0.0 

118.0 

1491.380 

0.0 

1.0 

4.55578E-06 

0.0 

173.0 

1490.780 

0.0 

1.0 

4.55578E-06 

0.0 

119.0 

1491.430 

0.0 

1.0 

4.55578E-06 

0.0 

174.0 

1490.760 

0.0 

1.0 

4.55578E-06 

0.0 

120.0 

1491.330 

0.0 

1.0 

4.55578E-06 

0.0 

175.0 

1490.720 

0.0 

1.0 

4.55578E-06 

0.0 

121.0 

1491.280 

0.0 

1.0 

4.55578E-06 

0.0 

176.0 

1490.710 

0.0 

1.0 

4.55578E-06 

0.0 

122.0 

1491.210 

0.0 

1.0 

4.55578E-06 

0.0 

177.0 

1490.740 

0.0 

1.0 

4.55578E-06 

0.0 

123.0 

1491.230 

0.0 

1.0 

4.55578E-06 

0.0 

178.0 

1490.720 

0.0 

1.0 

4.55578E-06 

0.0 

124.0 

1491.120 

0.0 

1.0 

4.55578E-06 

0.0 

179.0 

1490.760 

0.0 

1.0 

4.55578E-06 

0.0 

125.0 

1491.130 

0.0 

1.0 

4.55578E-06 

0.0 

180.0 

1490.780 

0.0 

1.0 

4.55578E-06 

0.0 

126.0 

1491.190 

0.0 

1.0 

4.55578E-06 

0.0 

181.0 

1490.800 

0.0 

1.0 

4.55578E-06 

0.0 

127.0 

1491.160 

0.0 

1.0 

4.55578E-06 

0.0 

182.0 

1490.810 

0.0 

1.0 

4.55578E-06 

0.0 

128.0 

1491.110 

0.0 

1.0 

4.55578E-06 

0.0 

183.0 

1490.810 

0.0 

1.0 

4.55578E-06 

0.0 

129.0 

1491.180 

0.0 

1.0 

4.55578E-06 

0.0 

184.0 

1490.740 

0.0 

1.0 

4.55578E-06 

0.0 

130.0 

1491.170 

0.0 

1.0 

4.55578E-06 

0.0 

185.0 

1490.730 

0.0 

1.0 

4.55578E-06 

0.0 

131.0 

1491.130 

0.0 

1.0 

4.55578E-06 

0.0 

186.0 

1490.690 

0.0 

1.0 

4.55578E-06 

0.0 

132.0 

1491.090 

0.0 

1.0 

4.55578E-06 

0.0 

187.0 

1490.640 

0.0 

1.0 

4.55578E-06 

0.0 

133.0 

1491.030 

0.0 

1.0 

4.55578E-06 

0.0 

188.0 

1490.610 

0.0 

1.0 

4.55578E-06 

0.0 

A-2 


1500  0.0  1028.0 

228.00  1881.00 
1028.00  3245.80 

>A’,  0. 

1028.00,  5200.00 
1490.0  5000.00 

15.0 

100  2.0  200.0  / 

48  101.375  189.5  / 


0.0 

2.06 
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0.0 
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0.0 

o 

o 
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APPENDIX  B 

SPECTRAL  ANALYSIS  CODE 


The  two  MATLAB*.m  files  shown  in  this  appendix  comprise  the  2-D  spectral  coherence  estima¬ 
tion  algorithm  used  to  generate  the  results  in  section  3.  The  function  make_index  is  first  called  to 
construct  the  INDX  matrix  used  to  map  between  rectangular  and  polar  coordinates.  Then  bispec  is 
called  to  generate  the  2-D  spectral  coherence  function. 

ftmction  [Indi.spec.size]  =  inake_iTidei(low_freq,high_freq,sr ,fft_3ize) 

7.  Function  ma)ce_index 

7.  [INDX,SPEC_SIZE]  =  MAKE_INDEX(LDW_FREq,HIGH_FREQ,SR,FFr_SIZE) 

Computes  the  mapping  matrix  INDX  which  is  used  in  the  estimation 
'L  of  spectral  coherence.  INDX  maps  points  on  a  rectangular  grid 

to  polar  coordinates. 

7. 

'!,  Variables 

X  L0W_FREQ:  Lowest  frequency  in  spectral  coherence  plane 

X  HIGH_FR£Q:  Highest  frequency  in  spectral  coherence  plane 

X  SR:  Sample  rate 

X  FFT.SIZE:  Length  of  FFT 

7. 

X  INDX:  Mapping  matrix 

X  SPEC_3IZE:  Used  to  determine  number  of  unique  slopes 

X  Author:  M.  Reuter  9-5-95 

X  Set  parameters  and  do  sanity  check 
if  low_freq  >=  high.freq 

dispC  ERROR:  Requested  low  frequency  is  greater  than  the  high  frequency') 
return 
end 

all_size  =  fft_size/2; 

low_index  =  round  (  (low_freq*2/sr)  *all_size)+l 
low_freq  =  3r*(low_index-l)/(2*all_size) 
high_index  =  round( (high_freq*2/3r) *all_size)+l 
high_freq  =  3r*(high_index-l)/(2*all_size) 
up_indei  =  roundC (high_freq*2/sr)*all_3ize)+l 

ii  UD  indejr->_all  -si  zp 

disp(’  Using’ ,sr/2, ’as  high  frequency’) 
up_indei  =  all_size 

end 

if  low_indei  <  1 

disp(’  Using  DC  as  low  frequency’) 
low_indei  =  1 

end 

m  =  up_ index  -  low.index  +  1; 

X  Allocate  memory 
X  =  zeros(m,m) ; 
polar_spec  =  zeros(m,2*m-l)  ; 
spec  =  zero3(m,m); 

Indi  =  zeros(2»m-l,m) ; 
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X  Generate  phasors  corresponding  to  each  coordinate 
Xr  =  kronC (low_indei-l : low_indei+m-2) .ones (m, 1) ) ; 

Xi  =  kxon(  (low_indei-l :  low_indei+iii-2)  ’  ,  ones  ( 1  ,m)  )  ; 

X  =  Xr  +  i*Xi; 

Ang  =  angle (X); 

spec.size  =  round(l+all_si2e-all_size/tan(Ang(m, 1) ) ) -1 ; 
clear  X  Xr  Xi 

%  Determine  indexes  corresponding  to  each  polar  coordinate 
for  kk=l : spec_size , 

angl  =  atan(all_size/ (all_size-kk+l) ) ; 

[Y,I]  =  min(abs(Ang  -  angl*ones(m,m) ) ) ; 

Indx(kk,l;m)  =  1(1 :m); 

end 

for  kk=l : spec_size-l , 

angl  =  atan((all_size-kk)/all_size) ; 

[Y,I]  =  min(abs(Ang  -  angl*ones (m,m) ) ) ; 

Indx(kk+spec_size ,  1  :m)  =  I(l:m); 

end 

high_f req , sr , L) 

7,  Function  bispec.indx 

7.  [SPEC]  =  BISPEC(INDX,SPEC_SIZE,CHAN1_SPEC,CHAN2_SPEC,L0W_FREQ, .  .  . 

7.  HIGH_FREQ,SR,L) 

7.  Computes  the  2D  spectral  coherence  function  of  two  time  series  by 

7.  employing  spectral  smoothing  using  a  rectangular  window. 

The  algorithm  first  maps  the  raw  spectrum  from  rectanguleir 
7.  coordinates  to  polcir  coordinates.  It  then  spectrally  smooths  the 

%  raw  spectrum  using  a  2D-convolution  with  the  rectangular  window. 

7  It  then  maps  the  smoothed  spectrum  back  to  rect^lngular  coordinates. 

7  INDX  and  SPEC_SIZE  must  have  been  generated  using  MAKE_INDEX. 

7 

7  Variables 

7  INDX:  Mapping  matrix 

*/  T-^T-*  .  tt__j  j..  ^  .  j. _ •  _  .  _ i- _  .n  _  n 
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low_index  =  round( (low_freq*2/sr)*all_size)+l 
loH_freq  =  sr*(loH_index-l)/(2*all_size) 
high_indei  =  roimd(  (high_freq*2/sr)  *all_size)  +  l 
high_freq  =  sr* (high_index-l) / (2*all_size) 
up.index  =  roundC (high_freq*2/sr)*all_size)+l 

if  up_index  >  all_size+l 

disp(’  Using’ ,sr/2 , 'as  high  frequency’) 
up_index  =  all.size 

end 

if  low_index  <  1 

disp(’  Using  DC  as  low  frequency’) 
low_index  =  1 

end 

m  =  up.index  -  low_index  +  1 ; 
y.L  =  roundCO .  020*all_size) 

Allocate  memory 
X  =  zeros(m,m) ; 
polar_spec  =  zeros(m,2*m-l)  ; 
spec  =  zeros(m,m); 

raw_spec  =  chanl_spec(loH_indei:up_index)*chan2_spec(lo»_indei:up_index) ’ ; 

%Hap  to  polar 
for  kk=l : 2*spec_size-l , 
for  11=1 :m 

polar_spec(ll ,kk)  =  raw.specCll ,Indx(kk , 11) )  ; 

end 

end 

7.  Spectral  smoothing  using  the  rectangular  window, 
polar.spec  =  conv2(  polax.spec  ,  ones(L,l),  ’same’  ); 

7ipolar_spec  =  conv2(  poleir_spec  ,  ones (L, 3),  ’same’  ); 

y.Map  to  rectangular 
for  kk=l : 2*spec_size-l , 
for  11=1 :m 

specCll ,Indx(kk , 11) )  =  polar_spec(ll,kk) ; 

end 

end 

y.Normalize  by  PSDs 

rpsd  =  diag(chcinl_spec(low_index;up_index)*chanl_spec(low_index;up_index)  ’)  ; 
Psdl  =  conv2(rpsd,  ones(L,l),  ’same’); 

rpsd  =  diag(chan2_spec(low_index:up_index)*chan2_spec(low_index:up_index) ’) ; 
Psd2  =  conv2(rpsd,  ones(L,l),  ’same’); 
spec  =  spec./(sqrt(Psdl*Psd2’)) ; 


B-3 


Figure  18. 
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